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Previous experiments have revealed that shock waves in thermally relaxing gases, such as ion- 
izing, dissociating and vibrationally excited gases, can become unstable. To date, the mechanism 
controlling this instability has not been resolved. Previous accounts of the D'yakov-Kontorovich in- 
stability, and Bethe-Zel'dovich-Thompson behaviour could not predict the experimentally observed 
instability. To address the mechanism controlling the instability, we study the propagation of shock 
waves in a simple two-dimensional dissipative hard disk molecular model. To account for the energy 
relaxation from translational degrees of freedom to higher modes within the shock wave structure, 
we allow inelastic collisions above an activation threshold. When the medium allows finite dissipa- 
tion, we find that the shock waves are unstable and form distinctive high density non-uniformities 
and convective rolls on their surface. Using analytical and numerical results for the shock Hugoniot, 
we show that both DK and BZT instabilities can be ruled out. Instead, the results suggest that the 
clustering instability of Goldhirsch and Zanetti in dissipative gases is the dominant mechanism. 
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Shock waves are well known to become unstable when 
travelling through an exothermic medium pQ. However, 
locally endothermic processes have also been shown to 
lead to shock instabilities, although their origin and phys- 
ical mechanism remain unclear. For example, sufficiently 
strong shock waves leading to ionization have been shown 
to develop instabilities [5]-[5] . Shock instabilities have also 
been observed in sufficiently strong shock waves when 
dissociation of the molecules is possible [S] . Perhaps the 
most puzzling experimental observations correspond to 
shock instabilities in carbon dioxide [SHI] and other larger 
molecules, in which neither ionization nor dissociation is 
expected. 

At present, an explanation for these observed instabilities 
is lacking. This is the focus of the present Letter. Current 
known accounts for shock instabilities in non-reactive 
media are the D'yakov-Kontorovich (DK) instability [5], 
which leads to undamped emission of sound from the 
shock and vortex-entropy waves, and the anomalous be- 
haviour in Bethe-Zel'dovich-Thompson fluids 9 which 
leads to shock splitting. Griffith et al. [5j attempted 
to rationalize their observations of shock instabilities in 
argon and carbon dioxide using the DK mechanism by 
quantitatively evaluating the DK instability criteria us- 
ing both experimental measurements of the shock Hugo- 
niot and equilibrium models. They however observed 
shock instabilities well in the region predicted to be sta- 
ble. Shock instabilities have also been attributed to hy- 
drodynamic instabilities in the low isentropic exponent 
gases in the presence of large shears, such as those due 
to large shock curvatures [7j. It is however difficult to 
reconcile this explanation with the instability observed in 
planar shock tube experiments. In spite of these difficul- 
ties to account for the experimental observations, shock 
instabilities have been reproduced in numerical experi- 
ments for ionizing shocks in argon by direct numerical 



simulation jTUJ [TT] . These simulations attribute the in- 
stabilities to the coupling of the nonlinear kinetics of the 
collisional-radiative model with wave propagation within 
the induction zone, akin to the mechanism of cellular for- 
mation in gaseous detonations pQ. Although the physi- 
cal mechanism remains unclear, these numerical findings 
suggest that the observed instabilities in previous exper- 
iments are not due to an experimental artefact, such as 
the contact surface Taylor instabilities prevalent in ex- 
periments. 

The present study attempts to elucidate the mechanism 
at play common to these previous experiments. To for- 
mulate our hypothesis, we first note that the common 
feature of the unstable configurations summarized above 
are that the compressed medium experiences a dissi- 
pative process during the relaxation to the equilibrium 
post-shock state. In strong shock waves, the first inter- 
nal modes to equilibrate within the shock structure are 
the translational and rotational degrees of freedom [T^] , 
within a few mean free paths. If the shock wave is suf- 
ficiently strong, vibrational modes can also be excited, 
but these take a significantly larger number of mean free 
times to equilibrate [12] . Likewise, if the temperature 
is sufficiently high to permit ionization and dissociation, 
the energy available in the translational and rotational 
modes is slowly patly transferred to these other internal 
modes through inelastic collisions. All these energy re- 
laxation processes correspond to an interval during which 
the translational and rotational kinetic energies of the 
molecules are lost through inelastic collisions. It is this 
inelasticity in the collision process that we wish to model. 
Indeed, Goldhirsch and Zanetti have shown that a clus- 
tering instability is possible in initially uniform gases un- 
dergoing inelastic collisions [13] . In the present study, 
we wish to test whether the same instability mechanism 
also operates in dissipative gases undergoing shock wave 
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relaxation to equilibrium. 

In order to capture the key relaxation processes involved, 
and retain a sufficiently simple model which will permit 
us to gain further analytical insight, we assume a two- 
dimensional medium composed of hard disks. We choose 
hard disks such that the simulations can be sufficiently 
simple to implement, test, analyze, visualize and repro- 
duce. Furthermore, we have previously validated an an- 
alytical formulation of the shock jump conditions in a 
hard disk gas [2]. This also permits us to analytically 
evaluate whether the DK and BZT mechanisms can ac- 
count for the instability. 

The model assumes that each binary collision is elastic, 
unless an activation threshold is reached. Activation is 
assumed to occur when the collision between two disks is 
sufficiently strong. This mimicks the excitation of higher 
degrees of freedom (rotation, vibration, dissociation, ion- 
ization, etc..) with increasing temperatures |12j . Quan- 
titatively, the collision between two disks is assumed to 
be elastic if the relative speed between two disks taken 
along the line of action is below a threshold u* , a classical 
activation formalism in chemical kinetics. For collisions 
with a higher amplitude, we assume an inelastic dissipa- 
tive collision modelled with a constant coefficient of resti- 
tution e < 1. The problem we study is a classical shock 
propagation problem, whereby the motion of a suddenly 
accelerated piston driven in a thermalized medium drives 
a strong shock wave. The piston is initially at rest and 
suddenly acquires a constant velocity u p . This model 
allows for the dissipation of the non-equilibrium energy 
accumulated within the shock structure which terminates 
once the collision amplitudes fall back below the activa- 
tion threshold. In this manner, the activation threshold 
also acts as a tunable parameter to control the equilib- 
rium temperature in the post shock media, see below. 
Note that the model assumed is also the standard model 
for granular gases [15] . and our results may find applica- 
tions in that field. 

The numerical experiments thus reconstruct the dynam- 
ics of the hard disks. These are calculated using the 
Event Driven Molecular Dynamics technique first intro- 
duced by Alder and Wainright 16 . We use the im- 
plementation of Poschel and Schwager [17], that we ex- 
tended to treat moving walls. The code was tested for 
non-dissipative media in our previous study |14j . The nu- 
merical experiments were performed using 30,000 disks, 
unless otherwise noted. The particles were initialized 
with equal speed and random directions. The system was 
let to thermalize and attain Maxwcll-Boltzmann statis- 
tics. Once thermalized, the piston started moving with 
constant speed. Below, all distances have been normal- 
ized by the radius of the disks and speeds by the initial 
speed of the disks. The initial packing factor of the disks 
was chosen to be rji = (V a /m)/v — 0.012, where V a /m 
is the specific volume (area) of the hard disk and v the 
specific volume; the initial gas is thus in the ideal gas 



regime [14] . 

The model was found to predict shock instabilities. Fig- 
ure [T] shows the structure of the compressed gas and its 
macroscopic averaged properties for a piston propagating 
with Up = 7, with e = 0.95 and an activation threshold of 
u* = 3. For reference, the shock Mach number is approx- 
imately 6.6. As can be clearly seen, the shock acquires 
a bumpy structure of approximately 10 — 15Ai spacing, 
where Ai is the mean-free-path evaluated at the initial 
state. Animation of the dynamics of the particles showed 
the evidence of macroscopic convection rolls. These were 
confirmed by coarse-graining the particle dynamics in or- 
der to visualize the streamlines. As can be seen in Figure 
[l] the large bumpy structure corresponds to vortex rolls, 
where the material accumulates in the large scale high 
density clusters. We have performed calculations with 
different domain sizes and number of particles, but this 
did not change the manifestation of the instability nor 
its characteristic spacing. 

Figure [T]d shows the coarse-grained density, temperature 
and pressure within the shock layer, obtained at the same 
evolution time as figure [T^i. These were obtained by en- 
semble averaging and coarse-graining in strips of width 
0.5Ai. Note that density and temperature were obtained 
directly from the particle data, while pressure was esti- 
mated using the Helfand equation of state for a dense 
hard disk medium [18 . 

The one-dimensional average shock structure resembles 
the characteristic structure of relaxing media: the trans- 
lational temperature has a peak within the shock and 
then decays while the non-equilibrium shock state dissi- 
pates its energy. For reference, Figure [T]d also shows the 
temperature evolution in a shock propagating through an 
elastic medium at the same initial conditions and shock 
Mach number (i.e., same mass flux across the shock). 
As can be seen, the peak in the dissipative calculation 
corresponds closely to the equilibrium translational tem- 
perature. Our model thus recovers the classical picture of 
relaxation phenomena inside shock waves (see Ref. [5] for 
example), where the first modes to relax are the transla- 
tional ones. The subsequent drop in temperature in the 
dissipative calculation is due to the inelastic collisions, 
which eventually become elastic again once the majority 
of the collisions become de-activated. 

Further insight into the shock transition structure and 
instability mechanism can be deduced by constructing 
the equilibrium shock Hugoniot, that is the locus in the 
pressure- volume plane of possible end states obtained af- 
ter shock compression. Figure [2] shows our numerical 
results obtained for different shock strengths for elastic 
collisions (e = 1) and inelastic collisions (e = 0.95). Also 
shown are the Hugoniot curves obtained analytically, us- 
ing the Helfand equation of state for a hard disk medium 
[14] , For the dissipative medium, the Hugoniot predic- 
tion captures very well the numerical results, as expected 
by our selection of the activation threshold to ensure an 
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FIG. 1. a) Particle distribution and coarse-grained stream- 
lines; b) coarse grained one-dimensional distribution of pres- 
sure, temperature and density for u p — 7, e = 0.95 and 
u* = 3; broken line is the shock temperature distribution 
for e=l and the same shock Mach number. 



approximately isothermal system. Also shown in Figure 
[2] is the evolution of the medium's state from initial to 
equilibrium state, for conditions illustrated in Figure [T] 
this transition line is termed the Rayleigh line. 

The shock Hugoniot and Rayleigh lines configurations 
obtained can immediately rule out any BZT behavior 
[9]. This would require a convex Hugoniot curve in the 
vicinity of the transition from elastic to inelastic behav- 
ior. The Hugoniot remains concave. Furthermore, BZT 
behavior would require the formation of a double discon- 
tinuity [9 , which again is not the case; the transition 
from initial to equilibrium state is achieved across a sin- 
gle discontinuity, since the Rayleigh line does not exhibit 
any convex behaviour. Since the slope of the line gives 
the mass flow rate across the shock wave, the constancy 
of the mass flow rate signifies a single discontinuity with 
a unique speed. Note however the slight concavity of the 
Rayleigh line, a signature of viscous stresses in the shock 




FIG. 2. The equilibrium states achieved across a system un- 
dergoing elastic (black symbols) and activated inelastic col- 
lisions (open symbols) and comparison with the ideal Hugo- 
niot and the isotherm computed with the Helfand equation of 
state. 



wave [9J. 

We can also determine whether DK instability is ex- 
pected or not based on the equilibrium shock Hugoniot 
data. The necessary condition for DK instability requires 
that the D'yakov parameter, h, be greater than a critical 
value h c , given by 



h = f 



dV 
dp 



> h, = 



H 



1 - a/|(i + Vi/v 2 ) 

1 - Mf (1 - Vi/V 2 ) 



(1) 



where the derivative is taken along the equilibrium shock 
Hugoniot and estimated at the equilibrium state 2, while 
j 2 = (p2 — Pi)/(Vi — V2) is the square of the mass flux 
density across the shock. For the dissipative hard disk 
gas considered, an isothermal equilibrium state approx- 
imates very well the Hugoniot. Using the isotherm to 
evaluate the derivative, we find that the DK criterion is 
never satisfied. The results are shown in Figure [3j Nei- 
ther the elastic nor the dissipative medium is expected 
to develop DK instabilities. We can thus conclude that 
neither BZT behavior nor the DK instability are com- 
patible with the results of the numerical experiments. 

To investigate whether the instability observed is com- 
patible with the clustering instability of homogeneous 
dissipative gases [T3], we conducted a separate simula- 
tion where all the gas is initialized with an average speed 
of u — 7 and random directions. This velocity was cho- 
sen to correspond to the piston speed investigated above. 
Once initiated, we did not include an activation thresh- 
old, and let the system cool via the homogeneous cooling 
state of granular dynamics (see [H]). We found that large 
scale density clusters appear on time scales only a few 
times longer than observed for the piston driven shock 
flows described above. While the time scales of the clus- 
tering formation in the homogeneous cooling state were 



FIG. 3. Propensity for D'yakov-Kontorovich instability in the 
ideal and in the dissipativc hard disk gas; h c — h < indicates 
instability. 



FIG. 4. Shocked layer morphology in a dissipative gas with 
a) e = 0.8 and b) e = 0.93. 



indeed expected to be longer than in the piston driven 
case, since the piston continuously energizes the shocked 
gas, the similarity in the time scales suggest that the 
same clustering instability is at hand. The large density 
clusters are formed by the inability of the gas within the 
cluster to leave owing to the local high rates of energy 
dissipation [13] , 

To gain further insight into the role of the rate of dis- 
sipation in the instability mechanism, we have further 
conducted numerical experiments by varying e. Varying 
e does not change the equilibrium final state, and does 
not affect the DK predictions. It only changed the rate 
at which energy is dissipated [TS] . Figure [4] shows the 
shock structure obtained with different values of e. We 
find that e controls both the relaxation length and the 
spacing of the bumps, which, to first approximation, was 
found to be approximately given by the relaxation length 
scale. This further supports our conclusion that the in- 
stability is associated with the rate of dissipation leading 
to cluster formation [15] . 

The numerical experiments conducted have thus shown 
that shock instabilities in relaxing media can be ac- 
counted for via the Goldhirsh and Zanetti clustering 
instability mechanism in dissipative gases, when the 
D'yakov-Kontorovich and Bethe-Zel'dovich-Thompson 
behavior can be ruled out. However, a quantitative pre- 
diction of the experimental observations in relaxing gases 
requires a much more realistic relaxation model, which we 
leave for further study. 
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